Linearized Dynamics Near Steady-State

We start from the linearized state-equations of the passively mode-locked laser ref1



In [1]:

    
import sympy as sym
sym.init_printing(use_latex='mathjax')

Pst, gst, EsatL, dqPdEP = sym.symbols("P_{st} g_{st} E_{satL} q'_P")
TR, tauL, etaP, Toc, PP0 = sym.symbols('T_R tau_L eta_P T_{oc} P_{P0}')
alpha, beta, gamma, epsilon = sym.symbols('alpha beta gamma epsilon')
w0 = sym.symbols('omega_0')

Pst_ = sym.Eq(Pst, -EsatL / tauL + etaP * PP0 / gst)

x = sym.MatrixSymbol('x', 2, 1)
xdot = sym.MatrixSymbol('xdot', 2, 1)
y = sym.MatrixSymbol('y', 1, 1)
u = sym.MatrixSymbol('u', 1, 1)

A = sym.Matrix([[-Pst * dqPdEP, Pst / TR],
                [-gst / EsatL, -1 / tauL - Pst / EsatL]])
B = sym.Matrix([[0], [etaP / EsatL]])
C = sym.Matrix([[Toc, 0]])

sym.Eq(xdot, A * x + B * u)









    Out[1]:





$$\dot{x} = \left[\begin{matrix}0\\\frac{\eta_{P}}{E_{{satL}}}\end{matrix}\right] u + \left[\begin{matrix}- P_{{st}} q'_{P} & \frac{P_{{st}}}{T_{R}}\\- \frac{g_{{st}}}{E_{{satL}}} & - \frac{1}{\tau_{L}} - \frac{P_{{st}}}{E_{{satL}}}\end{matrix}\right] x$$



In [2]:

    
sym.Eq(y, C * x)









    Out[2]:





$$y = \left[\begin{matrix}T_{{oc}} & 0\end{matrix}\right] x$$

where $\vec{x} = \left[ \delta P, \delta g \right]^T$, $\vec u = \left[\delta P_P \right]$, and $\vec y = \left[ \delta P_{out} \right]$.

Now, we rewrite this in new state variables $\vec x' = \left[ \delta \dot{P} / \omega_0, \delta P \right]$, with

$\omega_0 = \sqrt{ \frac{P_{st} g_{st}}{T_R E_{satL}} + P_{st} q'_P \left( 1 / \tau_L + P_{st} / E_{satL} \right)}$.

From the state equation above, we immediately see the transformation matrix $M$:



In [3]:

    
M = sym.Matrix([[-Pst * dqPdEP / w0, Pst / TR / w0],
                [1, 0]])
w0_ = sym.Eq(w0, sym.sqrt(-A[0, 1] * A[1, 0] + A[0, 0] * A[1, 1]))
xp = sym.MatrixSymbol("x'", 2, 1)
sym.simplify(sym.Eq(xp, M * x))









    Out[3]:





$$x' = \left[\begin{matrix}- \frac{P_{{st}} q'_{P}}{\omega_{0}} & \frac{P_{{st}}}{T_{R} \omega_{0}}\\1 & 0\end{matrix}\right] x$$

The transformed state equation is:

$ \begin{align*} \dot{x'} &= M A M^{-1} x' + M B u \\ y &= C \cdot x = C M^{-1} x' \end{align*}$



In [4]:

    
M = M.subs(w0_.lhs, w0_.rhs)

Bp = M * B
Ap = M * A * M**-1
sym.Eq(sym.MatrixSymbol("A'", 2, 2), sym.simplify(Ap))









    Out[4]:





$$A' = \left[\begin{matrix}- P_{{st}} q'_{P} - \frac{1}{\tau_{L}} - \frac{P_{{st}}}{E_{{satL}}} & - \sqrt{\frac{P_{{st}} q'_{P}}{\tau_{L}} + \frac{P_{{st}}^{2} q'_{P}}{E_{{satL}}} + \frac{P_{{st}} g_{{st}}}{E_{{satL}} T_{R}}}\\\sqrt{\frac{P_{{st}} q'_{P}}{\tau_{L}} + \frac{P_{{st}}^{2} q'_{P}}{E_{{satL}}} + \frac{P_{{st}} g_{{st}}}{E_{{satL}} T_{R}}} & 0\end{matrix}\right]$$



In [5]:

    
sym.Eq(sym.MatrixSymbol("B'", 2, 1), sym.simplify(Bp))









    Out[5]:





$$B' = \left[\begin{matrix}\frac{P_{{st}} \eta_{P}}{E_{{satL}} T_{R} \sqrt{\frac{P_{{st}}}{E_{{satL}} T_{R} \tau_{L}} \left(T_{R} q'_{P} \left(E_{{satL}} + P_{{st}} \tau_{L}\right) + g_{{st}} \tau_{L}\right)}}\\0\end{matrix}\right]$$

We recognize that the off-diagonal elements of $A$ are $\omega_0$ and $-\omega_0$. Furthermore, we define

$ \begin{align*} \zeta &= \frac{1}{2 \omega_0} \left(1 / \tau_L + P_{st} \cdot \left(q'_P + 1 / E_{satL} \right) \right) \\ \rho &= \frac{P_{st} \eta_P}{E_{satL} T_R \omega_0^2} \end{align*} $



In [6]:

    
zeta = sym.Symbol('zeta', real=True)
zeta_ = sym.Eq(zeta, (1 / tauL + Pst * (dqPdEP + 1 / EsatL)) / 2 / w0)
rho = sym.Symbol('rho', real=True)
rho_ = sym.Eq(rho, Pst * etaP / (EsatL * TR * w0) / w0)


Bpn = sym.Matrix(Bp)
Apn = sym.Matrix(Ap)

Bpn[0, 0] = Bp[0, 0] / (w0_.lhs / w0_.rhs) * rho_.lhs / rho_.rhs

Apn[0, 0] = Ap[0, 0] * zeta_.lhs / zeta_.rhs
Apn[0, 1] = Ap[0, 1] * w0_.lhs / w0_.rhs
Apn[1, 0] = Ap[1, 0] * w0_.lhs / w0_.rhs

sym.Eq(sym.MatrixSymbol("x'dot", 2, 1), sym.simplify(Apn * xp + Bpn * u))









    Out[6]:





$$\dot{x'} = \left[\begin{matrix}\omega_{0} \rho\\0\end{matrix}\right] u + \left[\begin{matrix}- 2 \omega_{0} \zeta & - \omega_{0}\\\omega_{0} & 0\end{matrix}\right] x'$$

Looking at the eigenvalues (see below), we see that the system

is stable if $\zeta > 0$
is unstable if $\zeta < 0$ (always)
has oscillatory modes for $-1 < \zeta < 1$ (radicand is negative)
is critically damped for $\zeta = 1$
is over-critically damped for $\zeta > 1$ (radicand is positive)

Of course, this analysis is only true if $\omega_0$ is real.



In [7]:

    
Apn.eigenvals()









    Out[7]:





$$\left \{ \omega_{0} \left(- \zeta - \sqrt{\left(\zeta - 1\right) \left(\zeta + 1\right)}\right) : 1, \quad \omega_{0} \left(- \zeta + \sqrt{\left(\zeta - 1\right) \left(\zeta + 1\right)}\right) : 1\right \}$$

Now, we apply a unit-step to the input of the system and determine where the output will stabilize. To do that, we have to solve the following equation:



In [8]:

    
eq = sym.Eq(sym.Matrix([[0], [0]]), Apn * xp + Bpn)
sym.simplify(eq)









    Out[8]:





$$\left[\begin{matrix}0\\0\end{matrix}\right] = \left[\begin{matrix}\omega_{0} \rho\\0\end{matrix}\right] + \left[\begin{matrix}- 2 \omega_{0} \zeta & - \omega_{0}\\\omega_{0} & 0\end{matrix}\right] x'$$

From the resulting $x'$ we obtain for the steady-state $y$ response:



In [9]:

    
sym.simplify(C * M**-1 * -Apn**-1 * Bpn)[0]









    Out[9]:





$$T_{{oc}} \rho$$

The DC-response to small pump-power changes is the slope efficiency.



In [10]:

    
sym.simplify((Toc * rho_.rhs).subs([(w0_.lhs, w0_.rhs),
                                    (Pst_.lhs, Pst_.rhs)]))









    Out[10]:





$$\frac{T_{{oc}} \eta_{P} g_{{st}}}{P_{{P0}} T_{R} \eta_{P} q'_{P} + g_{{st}}^{2}}$$

For $q'_P = 0$ (no saturable absorber), this simplifies to $\eta_P \cdot \frac{T_{oc}}{g_{st}}$, i.e. the slope efficiency depends only on the pump-absorption efficiency and the ratio of output coupling to total losses.